function z = respfct2Dexample(x)
% x contains currents, keff, and current eigenvalue
%   both j and keff are seeded by red-black


global order numg numel

%compute the response matrices
% 1 element case
if numel == 1
mm = [ 1 0 0 0 
       0 1 0 0
       0 0 1 0 
       0 0 0 1 ];                                % mm maps faces to faces
else
% 2 element case
mm = [ 1 0 0 0 0 0 0 0
       0 0 0 0 2 0 0 0
       0 0 1 0 0 0 0 0
       0 0 0 1 0 0 0 0
       0 2 0 0 0 0 0 0
       0 0 0 0 0 1 0 0
       0 0 0 0 0 0 1 0
       0 0 0 0 0 0 0 1];      
end

% tmp goes where reflection occurs
tmp  = diag( ones(1,order+1).*(-1).^(1+(1:order+1)) );
tmp  = kron( tmp, eye(numg) );
% tmp goes where transmission occurs
tmp2 = diag( ones(1,order+1).*(1).^(1+(1:order+1)) );
tmp2 = kron( tmp2, eye(numg) );
% M combines everything
M   = kron(mm==1, tmp )+kron(mm==2, tmp2);
M   = sparse(M);

% get the response functions
[R1 F1 A1 L1] = resp2dexample( x(end-1), order, 1  );
R1 = sparse(R1);
if (numel==1)
    R = R1; F = F1; A = A1; L = L1;
else
    % if two element
    [R2 F2 A2 L2] = resp2dexample( x(end-1), order, 2  );
    R2 = sparse(R2);
    R = sparse( [R1 zeros(size(R1)); zeros(size(R2)) R2] );
    F = [F1 F2];
    A = [A1 A2];
    L = [L1 L2];
end

j1 = M*x(1:end-2);

L=0;

    
%compute the nonlinear function, FF(x)*x = 0

FF = [ R*M-x(end)*eye(length(R(:,1)))  zeros(length(R(:,1)),1) 
       F*M                      -sum(A*M*x(1:end-2))-L   ];  
z1 = FF*(x(1:end-1)); 
z2 = 0.5 - 0.5*sum(x(1:end-2));

z  = [z1' z2']';


